! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2024-03-24: Initial creation.

module adv_2d_test_case_mod

  use adv_2d_common_mod

  implicit none

  real, parameter :: xmin = 0
  real, parameter :: xmax = 1
  real, parameter :: ymin = 0
  real, parameter :: ymax = 1
  real, parameter :: xrange = xmax - xmin
  real, parameter :: yrange = ymax - ymin
  real, allocatable :: x (:)              ! Cell center coordinates
  real, allocatable :: y (:)              ! Cell center coordinates
  real, allocatable :: xi(:)              ! Cell interface coordinates
  real, allocatable :: yi(:)              ! Cell interface coordinates
  real, allocatable :: u (:,:)            ! Velocity component along x axis
  real, allocatable :: v (:,:)            ! Velocity component along y axis
  real, allocatable :: ui(:,:)            ! Velocity component along x axis
  real, allocatable :: vi(:,:)            ! Velocity component along y axis
  real dx                                 ! Cell interval along x axis
  real dy                                 ! Cell interval along y axis
  real :: dt = 0                          ! Time step size
  integer :: nx = 100                     ! Cell number along x axis
  integer :: ny = 100                     ! Cell number along y axis
  integer nt                              ! Integration time step number
  real, parameter :: omega = 0.03         ! Rotation angular speed
  real, parameter :: x0 = 0.25            ! Initial coordinate
  real, parameter :: y0 = 0.5             ! Initial coordinate

contains

  subroutine adv_2d_test_case_init(case_name, ns, rho)

    character(*), intent(in) :: case_name
    integer, intent(in) :: ns
    real, intent(out) :: rho(1-ns:nx+ns,1-ns:ny+ns)

    integer i, j

    if (dt <= 0) then
      write(*, *) '[Error]: dt must be greater than 0!'
      stop 1
    end if
    nt = 2 * pi / omega / dt

    allocate(x (1-ns:nx+ns))
    allocate(y (1-ns:ny+ns))
    allocate(xi(1-ns:nx+ns))
    allocate(yi(1-ns:ny+ns))
    allocate(u (1-ns:nx+ns,1-ns:ny+ns))
    allocate(v (1-ns:nx+ns,1-ns:ny+ns))
    allocate(ui(1-ns:nx+ns,1-ns:ny+ns))
    allocate(vi(1-ns:nx+ns,1-ns:ny+ns))

    ! Set mesh grid coordinates.
    dx = xrange / nx
    do i = 1 - ns, nx + ns
      x(i) = (i - 1) * dx
    end do
    do i = 1 - ns, nx + ns
      xi(i) = x(i) + 0.5 * dx
    end do
    dy = yrange / ny
    do j = 1 - ns, ny + ns
      y(j) = (j - 1) * dy
    end do
    do j = 1 - ns, ny + ns
      yi(j) = y(j) + 0.5 * dy
    end do

    ! Set initial condition.
    do j = 1 - ns, ny + ns
      do i = 1 - ns, nx + ns
        u (i,j) = -omega * (y (j) - 0.5)
        v (i,j) =  omega * (x (i) - 0.5)
        ui(i,j) = -omega * (yi(j) - 0.5)
        vi(i,j) =  omega * (xi(i) - 0.5)
      end do
    end do

    select case (case_name)
    case ('slotted_cylinder')
      do j = 1, ny
        do i = 1, nx
          if (abs(x(i) - x0) >= 0.02 .and. (x(i) - x0)**2 + (y(j) - y0)**2 < 0.01) then
            rho(i,j) = 1.0
          else if (abs(x(i) - x0) < 0.02 .and. (x(i) - x0)**2 + (y(j) - y0)**2 < 0.01 .and. y(j) >= 0.55) then
            rho(i,j) = 1.0
          else
            rho(i,j) = 0.0
          end if
        end do
      end do
    case default
      write(*, *) '[Error]: Invalid case_name ', trim(case_name), '!'
      stop 1
    end select

    call apply_bc(ns, nx, ny, rho)

  end subroutine adv_2d_test_case_init

  subroutine adv_2d_test_case_final()

    if (allocated(x )) deallocate(x )
    if (allocated(y )) deallocate(y )
    if (allocated(xi)) deallocate(xi)
    if (allocated(yi)) deallocate(yi)
    if (allocated(u )) deallocate(u )
    if (allocated(v )) deallocate(v )
    if (allocated(ui)) deallocate(ui)
    if (allocated(vi)) deallocate(vi)

  end subroutine adv_2d_test_case_final

end module adv_2d_test_case_mod